bcftools使用笔记 |
您所在的位置:网站首页 › vcf文件怎么打开linux 命令 › bcftools使用笔记 |
![]() 什么是vcf?全称是variant call format,是用来记录SNP,INDEL和SV变异信息的文本。在GATK软件中得到最好的支持,当然SAMtools得到的结果也是VCF格式,和GATK的VCF格式有点差别。 1.查看vcf文件的头部信息bcftools view -h shrimp.vcf.gz 2.提取指定个体产生新vcf文件如果只想分析vcf文件中部分样本,可以利用view命令中的—S选项过滤个体。 bcftools view -S selectedinds.txt shrimp.vcf.gz -Ov > shrimp104.vcf -S选项后边跟一个文本文件,每一行为保留个体的ID编号。如果样品少,也可以 在-S 后边直接跟样品的ID号。>为重定向符号,表示把保留的个体信息存到shrimp.vcf文件中;-Ov表示输出未经压缩的vcf文件。 3.压缩vcf为gz压缩文件bcftools一般要求读入文件为gz格式。而且gz格式的文件,压缩幅度非常大,便于传输。因此可以利用view命令把vcf文件压缩为gz格式。 bcftools view shrimp.vcf -Oz -o shrimp.vcf.gz --threads 8 -Oz表示输出格式为压缩文件gz格式;-o后边跟压缩文件名字;–threads 8 表示通过8个线程并行压缩shrimp.vcf文件。 4.提取等位基因和基因型信息参考这个链接Extracting allele, Genotype from VCF file. 用到的命令为“query”,说明书中该命令的功能为“transform VCF/BCF into user-defined formats”,即把VCF/BCF文件转换为用户定义格式。 bcftools query -f '%CHROM %ID %POS %REF %ALT [ %TGT]\n' query.vcf.gz -o query.extract.txt 上述命令提取vcf文件中染色体、基因型等信息,输出为空格分隔的文本文件。 其中: %CHROM 染色体列%ID 变异位点名称%POS 变异位点位置%REF 参考等位基因%ALT 变异等位基因%TGT 字符格式如A/G的基因型;%GT为0/1格式的基因型也可以输出为英文逗号分隔的csv文件。若想用tab分隔,符号是bcftools query -f ‘%CHROM,%ID,%POS,%REF,%ALT[,%TGT]’ query.vcf.gz -o query.extract.csv` 5.变异位点的基本统计分析stats命令用于统计VCF文件的基本信息,比如突变位点的总数,不同类型突变位点的个数等。用法如下: bcftools stats view.vcf > view.stats 然后会生成一个名为view.stats的文本文件。 然后调用一下命令,进行可视化输出: plot-vcfstats view.stats -p output 多个结果文件保存在output文件夹下。其中summary.pdf文件,包括了汇总分析的结果。主要包括: 转换(Transition:嘌呤(purine)转换A(腺嘌呤)-G(鸟嘌呤);嘧啶(pyrimidine)转换C(胞嘧啶)-T(胸腺嘧啶))和颠换(Transversion:A-C;A-T;G-C;G-T)类型分布统计。尽管颠换类型要多于转换,但在基因组上,转换的实际发生频率要高于颠换类型。 一般转换与颠换的发生比例为2:1。需要注意,plot-vcfstats脚本在bcftools安装目录misc文件夹下,这是一个perl脚本,会调用python3绘图模块Matplotlib。如果没有安装该模块,可以通过pip3命令pip3 install -U matplotlib进行安装。如果没有pip3,通过sudo apt-get install python3-pip命令安装。 plot-vcfstats生成pdf文件还需要pdf-latex,如果系统没有安装latex,通过sudo apt-get install texlive-full进行安装。 上述命令在ubuntu 18.04系统中执行。 6. 替换染色体名称需要用到annotate命令中的–rename-chrs参数。 命令形式: bcftools annotate --rename-chrs NewChrName.txt old.xxx.vcf.gz -Oz -o new.xxx.gz --threads n 其中NewChrName.txt文件中包括了旧和新染色体名称的对应关系。--threads可以设置多线程加快新vcf文件的生成速度。 可以通过bioawk命令来提取旧染色体名字。 bioawk -c vcf ' { print $chrom } ' old.vcf.gz | sort -n | uniq > ChrName.txt |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |